Anomalous lateral diffusion in a viscous membrane surrounded by 
viscoelastic media 
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Abstract - We investigate the lateral dynamics in a purely viscous lipid membrane surrounded 
by viscoelastic media such as polymeric solutions. We first obtain the generalized frequency- 
dependent mobility tensor and focus on the case when the solvent is sandwiched by hard walls. 
Due to the viscoelasticity of the solvent, the mean square displacement of a disk embedded in the 
membrane exhibits an anomalous diffusion. An useful relation which connects the mean square 
displacement and the solvent modulus is provided. We also calculate the cross-correlation of the 
particle displacements which can be applied for two-particle tracking experiments. 
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Introduction. — Biomembranes are thin two- 
dimensional (2D) fluids which separate inner and outer en- 
vironments of organelles in cells. The fluidity of biomem- 
branes is guaranteed mainly due to the lipid molecules 
which are in the liquid crystalline state at physiological 
temperatures. Proteins and other molecules embedded in 
biomembranes undergo lateral diffusion which plays im- 
portant roles for biological functions [1]. It should be 
noted, however, that biomembranes are not isolated 2D 
systems, but are coupled to the surrounding polar solvent 
such as water. Indeed the presence of water is essential for 
amphiphilic lipid molecules to spontaneously form bilayers 
by self-assembly. 

In the last few decades, it has been recognized that 
the outer solvent has a significant effect on the mem- 
brane dynamics which takes place in 2D. SafFman and 
Delbriick considered the Brownian motion of a 2D disk 
confined in a fluid membrane [2l[3]- In their hydrody- 
namic model, the transfer of membrane momentum to the 
outer three-dimensional (3D) fluid was taken into account 
through the boundary conditions at the membrane sur- 
faces. The translational diffusion coefficient in the weak 
coupling (small disk) limit was shown to exhibit a loga- 
rithmic dependence on the disk size. This dependence has 
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been repeatedly tested for various lipid molecules and pro- 
teins [3]. In the strong coupling (large disk) limit [5], on 
the other hand, the diffusion coefficient is inversely pro- 
portional to the disk size showing the analogy to the 3D 
Stokes-Einstein relation [B]. Such a 3D-like behavior in 
2D membrane is caused by the back flow effect mediated 
by the bulk solvent. 

In this Letter, we discuss the dynamics and responses of 
membranes when their surrounding solvent is viscoelastic 
rather than purely viscous. This is a common situation in 
all eukaryotic cells whose cytoplasm is a soup of proteins 
and organelles, including a thick sub-membrane layer of 
actin-meshwork forming a part of the cell cytoskeleton [T] . 
The extra-cellular fluid can also be viscoelastic because it 
is filled with extracellular matrix or hyaluronic acid gel. 
In addition to the surrounding solvents, lipid membranes 
themselves can be viscoelastic [7] . Although it turned out 
that membranes are purely viscous in the latest report [7] , 
their experimental technique using particle tracking mi- 
crorheology provides us with a new clue to investigate the 
dynamical responses of lipid bilayers coupled to the sur- 
rounding environments under controlled conditions. Re- 
cently, viscoelasticity of phospholipid Langmuir monolay- 
ers in a liquid- condensed phase was measured using active 
microrheology [8]. Being motivated by these works, we 
discuss the mean square displacement (MSD) of a circular 
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disk embedded in a 2D sheet by taking into account the 
viscoelasticity of the surrounding media. This quantity is 
experimentahy measurable in single-particle tracking mi- 
crorheology [Hllin]- We further calculate the two-particle 
MSD which is useful for two-point microrheology experi- 
ments [TT]. For both cases, wc show that the viscoelastic- 
ity of the surrounding media leads to an anomalous diffu- 
sion in the 2D viscous membrane. 

Recently, Granek discussed the dynamics of an undu- 
lating bilayer membrane surrounded by viscoclastic me- 
dia [12]. He calculated the frequency-dependent trans- 
verse (out-of-plane) MSD of a membrane segment and the 
linear response to external forces. In our theory, we treat 
the membrane as an infinitely large flat sheet with a 2D 
viscosity rj, and consider its lateral dynamics. A simi- 
lar problem was considered in refs. [131114] in which the 
authors have taken into account the viscoelasticity of the 
membrane itself because they were originally motivated by 
the earlier experiment of ref. [7]. A more general theory 
for the dynamics of viscoelastic membranes was given by 
Levine and MacKintosh [TS]. In these works, however, the 
surrounding solvent is assumed to be purely viscous. Sim- 
ilar to Granek's work, the main purpose of our work is to 
point out the importance of the viscoelasticity of the bulk 
solvent on the membrane lateral dynamics such as diffu- 
sion or linear viscoclastic response. To make this point 
clear enough, we intentionally treat the membrane as a 
purely viscous 2D fluid. 

Hydrodynamic model. — We first establish the gov- 
erning hydrodynamic equations for our model. As shown 
in fig. [TJ the fluid membrane, fixed in the a;y-plane at 
2: = 0, is embedded in a bulk solvent that is further 
bounded by hard walls at z = ±/i. Let v(r, t) be the 
2D velocity of the membrane fluid at position r ~ (x, y) 
and at time t. We assume that the membrane is incom- 
pressible; 

V-v = 0. (1) 

We also work in the low-Reynolds number regime so that 
the inertial effects can be neglected. Then the Stokes equa- 
tion for the fluid membrane is given by 



dv 2 ^ 

p— 77V V - Vp 



F, 



(2) 



where p is the membrane 2D density, rj the membrane 
2D constant viscosity, p the in-plane pressure, fg the force 
due to the solvent described later, and F is any other force 
acting on the membrane. Notice that V stands for the 2D 
differential operator. 

Next we give the equations for the surrounding vis- 
coelastic solvent. The upper {z > 0) and the lower (z < 0) 
regions of the solvent are denoted by "+" and "— " , respec- 
tively. The velocities and pressures in these regions are 
written as v^(r,z,t) and p^{r,z,t), respectively. Similar 
to the fluid membrane, the solvent is also assumed to be 
incompressible ; 
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Fig. 1: Schematic picture showing a planar viscous membrane 
at 2 = with 2D constant viscosity t]. It is sandwiched by a 
viscoelastic solvent of 3D frequency-dependent viscosity 77s [t<j]. 
Two impenetrable walls are located at 2 = ±/i bounding the 
solvent. 



The Stokes equation for the viscoclastic solvent is written 
as 

P^^T- = V ■ o-^ 



dt 



(4) 



where is the solvent 3D density (assumed to be the same 
for both solvents) , and V indicates the 3D differential op- 
erator. Within the linear viscoelasticity approximation, 
the stress tensor cr^ in the above equation is given by the 
following constitutive relation [16] 



(5) 



0. 



(3) 



where ris{t) is the time-dependent solvent viscosity (as- 
sumed to be the same for both solvents). The rate-of- 
strain tensor is given by = [Vv^ -t- (Vv^)'^]/2, where 
the superscript "T" denotes the transpose. The surround- 
ing viscoelastic solvent exerts force fg on the membrane 
as considered in eq. It is given by the projection of 
((T+ — 0-^)2=0 ■ 62 on the xy-plane, where is the unit 
vector along the z-axis. We employ the stick boundary 
conditions at any time, i.e., the equalities of the veloci- 
ties, at z = and z = ±/i. When the solvent is purely vis- 
cous, this model reduces to that considered by the present 
authors [17] . 

It is convenient to perform the 2D Fourier transform 
in space and the Fourier-Laplace (or one-sided Fourier) 
transform in time for any function /(r,i) as dcflncd by 
/[k, w] = d^r dt f{r, t) exp[-i(k • r + wt)], where 
k ~ {kx,ky) is the 2D wavevector and w the angular fre- 
quency. Assuming that fluids are at rest at t = 0, we cal- 
culate is and obtain the membrane velocity as v[k, w] = 
G[k, oj] •F[k, w], where G[k, cj] is the frequency-dependent 
mobility tensor. Following the similar procedure described 



p-2 



in refs. [T31II1], we obtain 
Gaf} [k, w] 



rjk'^ + 2ris[i-o]k' coth{k'h) + iujp' 



(6) 



1/2 

where k' = k {l + iujps/ris[uj]k'^) , k = |k| and a, /3 = 
X, y. In the limit of ft. oo, which we call as the "free 
membrane case" , the above mobility tensor reduces to 



Sal3 ~ kakfj/P 



Tjk'^ + 2i-^s[^]kyJ\ + ia;ps/?7s[w]fc2 + iujp 



(7) 

which was given in ref. [13]. The opposite limit of ft. — > 
is called as the "confined membrane case" and the corre- 
sponding mobility tensor becomes 



5al3 - kgk/s/k'^ 



(8) 



where we have introduced K[a;] whose square is defined as 



rjh 



lujp 



(9) 



In the limit of w 0, eq. ^ reduces to the static mo- 
bility tensor obtained in ref. [17] . The free membrane case 
was originally considered by Saffman and Delbriick [2l|3], 
and the relevant length scale is the Saffman-Delbriick 
length ri/r]s beyond which the membrane feels the pres- 
ence of the outer solvent. On the other hand, the confined 
membrane case corresponds to that of a supported mem- 
brane close to the substrate as considered by Evans and 
Sackmann [18] and later generalized by us [TS]. Here the 
corresponding hydrodynamic screening length is set by the 
geometric mean of the Saffman-Delbriick length jy/jya and 
the distance ft between the membrane and the wall, i.e., 
[20]. Notice that the presence of the second wall 
only doubles the screening effect. In the following, we shall 
mainly consider the confined membrane case which allows 
us to treat most of the calculations analytically. This is 
mainly because k[ijj] does not depend on k. However, it 
should be noted that the limiting expression of eq. ^ 
gives a reasonable approximation to the full expression of 
eq. ([6]) even for ft ~ r]/r]s at least for a large enough moving 
object [201121]. 

Concerning the viscoelasticity of the surrounding sol- 
vent, we assume that its complex modulus obeys the 
power-law behavior such that Gs[w] = Go(ja;)" with 
a < 1, as generally argued by Granek [12]. This be- 
havior is commonly observed for various polymeric so- 
lutions at high frequencies. Examples are a = 1/2 and 
a = 2/3 for Rouse and Zimm dynamics, respectively [22) . 
a = 3/4 for semi-dilute solutions of semi- flexible poly- 
mers such as actin filaments |23j . From the viewpoint of 
particle-tracking microrheology experiment [24 1 l25 j . it is 
more convenient to work in the Laplace domain defined 
by /(s) = /q°° di /(i)e^''*. Since we have been working in 
the Fourier-Laplace domain, it is straightforward to con- 
vert to the Laplace domain by substituting s = itj. From 



Gs(s) = Gqs", the Laplace transform of the solvent vis- 
cosity behaves as 77s (s) = Gs(s)/s = Gqs"^^. Then cq. ^ 
simply becomes 



2Gos° 



rjh 



ps 



2Gos"~i 
77ft 



(10) 



where we have dropped the inertial term at the end. This 
approximation is justified when Gos""^/pft ^ 1, which 
is always valid for s — > corresponding to the long-time 
behavior. 

Single-particle tracking. — For the description of 
the Brownian motion of a circular disk confined in a mem- 
brane, we basically follow the formulation in ref. [26] . Let 
a and m be the radius and the mass of the disk, respec- 
tively. The effective generalized Langevin equation is writ- 
ten as [57] 



dt' \{t-t')U{t') + R{t), (11) 



where U{t) is the velocity of the disk, and X{t) is the time- 
dependent drag coefficient given below. The random force 
R{t) satisfies the fiuctuation dissipation theorem (FDT) 
when averaged over the ensemble of molecular motions; 
i.e., {R{t)) = and {R{to)R{to + t)) = k^TXit), where 
fee is the Boltzmann constant and T the temperature. In 
ref. j26j , it was shown that the renormalized mass is given 
by m* = m + irpa^ which takes into account the addi- 
tional inertia due to the drag from the fiuid membrane. 
Furthermore, the Laplace transform of X{t) is calculated 
to be 



A(.) 



2Trfjs{s)a'^ 
h 



47r?7 



k(s)aKi [k{s)a\ 
Ko[k{s)a] 



(12) 



where Kq[x] and A'i[x] are modified Bessel functions of 
the second kind, order zero and one, respectively. 

Following the standard procedure to obtain the 
MSD [27], one can relate its Laplace transform and A(s) 
through 



(Af^(.)) 



s2A(s)' 



(13) 



where the number of degrees of freedom tracked in the 
MSD is chosen to be two. For general t, application of 
the inverse Laplace transform provides us with the time- 
dependent MSD; 



1 



C+200 



, 4fcBT 
ds — = e . 

s2A(s) 



(14) 



In order to demonstrate how to use the above relations, 
we consider here the limit of a large disk size a — > 00 
{k{s)a ^ 1) so that the drag coefficient in eq. ([12]) is dom- 
inated by the first term, i.e., A(s) « 27Tfjs{s)a'^ / h. Then 
the Laplace transformed MSD in eq. ([T^ becomes 



2kBTh 



2kBTh 



(15) 
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This is the equation which relates the observed 2D MSD 
to the modulus of the surrounding bulk solvent (rather 
than the membrane). In other words, we can extract the 
solvent 3D information by using the 2D information due 
to the motion of a disk in the membrane. Once Gs{s) is 
obtained from the experiment, the frequency dependence 
of the storage and the loss moduli can be deduced by 
identifying Gs[lo] = G',[lo] + iG"[w] = Gs{s = iu). No- 
tice that these two representations are equivalent because 
Gi.[uj] and G'"[a;] are related by the Kramers-Kronig rela- 
tion [MIESII^. 

On the other hand, suppose the modulus of the bulk sol- 
vent is known a priori to behave as Gs{s) = Gqs" from in- 
dependent experiments, the MSD can be simply obtained 
from eqs. and ([T51) as 



2kBTh 



2^GoT[l 



(16) 



where T[x] is the gamma function. This calculation shows 
that the viscoelasticity of the solvent results in a subd- 
iffusive time dependence of the MSD. Since a < 1, the 
viscoelasticity slows down the normal diffusion process. 
This is the main result of this Letter. Compared with the 
3D case [9l[T0], the above expression is unique because it 
is proportional to h/a?. This 1/a^-dependence arises from 
the mass conservation in 2D rather than the momentum 
conservation [17]. 

In the limit of a small disk size a — > {'k{s)a <^ 1), 
the situation is more complicated. In this case, the drag 
coefRcient asymptotically behaves as 



A(s) w A-K-q 



In 



7 



(17) 



where 7 = 0.5772 
for k(s), a similar calculation yields 



k{s)a ^ 

is Euler's constant. By using eq. (HU]) 
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In 



f2r]ht° 



-I- (a — 3)7 — a + 1 



V Goa2 

(18) 

Since a < 1, this MSD grows like tln(l/i). Such a loga- 
rithmic correction leads to a time-dependent diffusivity. It 
should be noted, however, that this asymptotic expression 
is valid only when Goo^t^^" /rjh ^ 1. We also mention 
that all the above expressions recover our previous results 
when the solvent is purely viscous, i.e., a ^ 1 [19]. 

In fig. [2I we plot the result of numerical inverse 
Laplace transform of eq. ([T4]) with the full drag co- 
efficient eq. (fT2]) when a = 1/2, and compared 
it with the asymptotic expressions; eqs. (|T51) and 
(fT5|) . Here the dimensionlcss MSD and time are 
given by (27rr?/fcBT)(?7/i/2Goa2)i/("-i) (Ar2(t)) and t = 
{r}h/2Goa'^y^^°'~^h, respectively. As mentioned above, 
eq. (fT8)) is valid only for short time {i <^ 1), while 
eq. (fT6]) describes the asymptotic long-time behavior {t ^ 
16i/(i-")). The intermediate time region is described by 
an apparent power-law which will be explained 

elsewhere. 




Fig. 2: Scaled MSD as a function of scaled time i when a — 1/2. 
These dimensionless quantities are defined in the text. The 
dotted and dashed lines correspond to the asymptotic expres- 
sions given by eqs. H16|) and (I18|) . respectively. 



Two-particle tracking. — So far, we have discussed 
the motion of a single disk of finite radius a. As discussed 
in ref . [25] , there are several advantages to perform multi- 
particle microrheology. For example, long-time convective 
drift can be automatically subtracted in this method so 
that measurements of probe self-diffusivities become pos- 
sible over longer times. Multi-particle techniques can be 
also used to investigate heterogeneous materials. Here we 
discuss the effects of the solvent viscosity on the cross- 
correlation of two distinct probe positions, namely, two- 
point microrheology [11] . We show below that the distance 
between the two points essentially corresponds to the size 
of a disk in single-particle microrheology. 

Consider a pair of point particles embedded in the mem- 
brane undergoing Brownian motion separated by a 2D vec- 
tor r. The quantity of interest is the cross-correlation of 
the particle displacements {Ar^{0)Ar'p{t))r, where Ar^ is 
the displacement of the particle i (= 1,2) along the axis a 
(= X, y). We also define the x-axis to be along the line con- 
necting the two particles, i.e., r = re^- According to the 
FDT, this correlation function is related to the coupling 
mobility Ma[j{t) in the Laplace domain as |25) 

{ArlArl{s))r = -^M^p[r,s), (19) 

for sufficiently large r. The above equation is the analog 
of eq. (|13p for the two-particle tracking. Since M^y = by 
symmetry, we consider the longitudinal coupling mobility 
A/l = Mxx and the transverse one Mt = Myy. Note that 
the coupling mobility is not the inverse of the coupling 
resistance in the multi-particle case |25j . By utilizing the 
results in refs. [T7l[28l[29| . one can obtain these mobilities 
analytically. 
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First the Laplace transform of the longitudinal coupling 
mobility turns out to be 



1 



(20) 



In the limit of a large distance r — > oo (k(s)r ^ 1), the 
above expression can be approximated as 



ML(r,s) 



1 



1 



27rry {k{s)ry 



(21) 



As in the calculation of MSD of a single disk, we use 
eq. (|10p for k{s) and perform the inverse Laplace trans- 
form of eq. (jl9p . Then we obtain 



k^Th 



27rG'or[l + a] 



(22) 



The subdiffusive dependence on time and the 
dependence on distance is analogous to eq. ([TS]) , implying 
that r corresponds to a. 

In the limit of a small distance r — ^ {k{s)r ^ 1), on 
the other hand, eq. (|20)) asymptotically behaves as 



ML(r,s) 



1 



In 



K(s)r 



■7 + 



1 



(23) 



Following the same process as above, the cross-correlation 
function asymptotically behaves as 



kBT^ 



In 



f2r]ht° 



V Gor- 



(24) 



(a - 3)7 - a + 2 



which is valid for Gor^t^'" /-qh < 1. Notice that eq. ([21 
is also analogous to eq. (fT8|). 

Next the transverse coupling mobility is given by 



^T(r,s) = ^ 



i^o[K(s)r] + 



1 



Ki [K{s)r] 

fi{s)r (k(s)ry 



(25) 

The large distance limit (r — > 00) and the small distance 
limit (r — >■ 0) of this expression are 



and 



Mrir, s) 



MT(r,s) « — 



27r?7 (K(s)r)2 ' 





1" 











(26) 



(27) 



respectively. Since these forms differ from eqs. (PT|) and 
only by a sign, we do not repeat here the same cal- 
culations. However, as far as the time dependence of 
{Ary{0)Ary{t))r is concerned, it is essentially given by 
eqs. (gl]) and 



Discussion. — In this Letter, we have discussed the 
dynamics in a purely viscous lipid membrane surrounded 
by viscoelastic solvents such as polymeric solutions. Us- 
ing the generalized frequency-dependent mobility tensor 
for the confined membrane case, we calculated the MSD 
of a disk embedded in the membrane and obtained some 
asymptotic expressions. The obtained MSD exhibits an 
anomalous diffusion reflecting the viscoelastic property of 
the bulk solvent. For single-particle microrheology experi- 
ments, we presented an useful relation which connects the 
MSD and the solvent modulus in the Laplace domain when 
the size of the disk is large enough. We also obtained the 
cross-correlation of the particle displacements which can 
be used for two-particle tracking experiments. Our theory 
can be applied not only for lipid membranes but also for 
Langmuir monolayers. 

It is worthwhile to point out the implicit assumptions 
which are used in the present theory |25] . First we 
have assumed that the system obeys the FDT which al- 
lows us to relate the thermal fluctuations of the probe 
disk directly to the time-dependent drag coefficient (see 
eq. ([13])) or the coupling mobility (see eq. ((HI)). How- 
ever, in non-equilibrium situations with membrane trans- 
port proteins, for instance, the FDT can be violated and 
the deviation from the present theory may arise. Recently, 
it was reported that an in vitro model system consist- 
ing of a cross-linked actin network with embedded force- 
generating myosin motors strongly violates FDT [5D] . The 
second crucial assumption is that the generalized drag co- 
efficient or the coupling mobility are given precisely by 
their Newtonian analogs, but with the Newtonian viscos- 
ity r]s replaced by the frequency-dependent complex vis- 
cosity 77s [w] at all frequencies (see eq. (O). This is not 
an obvious assumption, but it has proved to be quite suc- 
cessful for any probe motion in a simple linear viscoelastic 
material under various conditions at least in 3D |25| . 

Some caution is required when applying our theory to 
experiments if the viscoelastic solvent is, for instance, a 
semi-dilute polymer solution. For single- particle tracking, 
our continuum approach is valid for inclusion sizes much 
larger than the mesh size of the network. This can be 
relevant such as for a micron-size membrane domain and 
a network correlation length of tens of nanometers. How- 
ever our theory may not be applicable for small membrane 
proteins which feel the solvent as purely viscous. For two- 
particle tracking, the distance between the two point par- 
ticles (which can be small membrane proteins) should be 
larger than the network mesh size. In order to obtain the 
full time behavior of the particle motion, one should use a 
viscoelastic modulus that is dependent both on wavevector 
and frequency. On the other hand, any simple fluid can 
be effectively viewed as viscoelastic at very high molecular 
frequencies [31] . Our model can be also used to study such 
a high frequency dynamics by using a Maxwell model for 
the solvent. 

In the present work, we have mainly discussed the con- 
fined membrane case in order to obtain analytical expres- 
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sions. Unfortunately, a single analytic expression of the 
drag coefficient for the whole range of the disk size is not 
known for the free membrane case. However, in the limit 
of a large disk size a — >■ cx), the asymptotic expression was 
obtained by Hughes et al. as Aftoc ~ Wrjsa which depends 
only on the solvent viscosity and is proportional to the 
size [5] . Assuming the replacement of the solvent viscosity 
with the time-dependent one as discussed above, the MSD 
of the disk becomes 



(Af2(s))ft.cc = 



(28) 



When the modulus of the bulk solvent behaves as Gs(s) — 
Gqs" as before, the MSD in the free membrane case is 
given by 



4aGol [1 + a\ 

which is again a subdifFusive behavior. In the opposite 
limit of a ^ 0, the Laplace transformed drag coefficient 
can be written as [HIS] 



Afrcc(s) ~ 47r?7 



It then follows that 

(Ar2(t))f,ee = 

TTT] 



In 



In 



Go a 



(30) 



(31) 



for Goat^~" /rj <^ 1. Although the above argument is 
not rigorous, it essentially captures the anomalous diffu- 
sion behavior in the free membrane case. For the two- 
particle tracking case, however, the corresponding cross- 
correlation function can be obtained without any ambi- 
guity since both the longitudinal and transverse coupling 
mobilities were analytically obtained for the free mem- 
brane case [29] . 

Currently we are extending our model to the case when 
the two solvents have asymmetric viscoelastic properties 
or when the disk is actively driven. Additional effects 
such as the finite curvature of vesicles or the out-of-plane 
deformation of the membrane should also be taken into 
account. 
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